Landsat SAV

Code

import numpy as np, pandas as pd, geopandas as gpd, sklearn.metrics as metrics
import geemap, ee , rasterio, sankee, warnings, pickle
from matplotlib.lines import Line2D
from rasterio import features
from geopandas import GeoDataFrame
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
from matplotlib.colors import ListedColormap, LinearSegmentedColormap
from rasterio.mask import mask
from osgeo import gdal, osr, gdalconst
import os, glob, subprocess, json, fiona
import numpy.ma as ma
import rasterio, pickle
import matplotlib as mpl
import matplotlib.pyplot as plt, plotly.express as px
from sklearn import ensemble
from sklearn.cluster import DBSCAN, AgglomerativeClustering
from matplotlib_scalebar.scalebar import ScaleBar
from scipy import stats
import seaborn as sns
import statsmodels
import glob
from glob import glob
import rasterio.plot as rplot
from scipy.cluster.hierarchy import dendrogram, linkage
from shapely.geometry import Polygon, Point, box
from sklearn.metrics import confusion_matrix, cohen_kappa_score
from rasterio.plot import show
import rioxarray as rxr
from rio_cogeo.profiles import cog_profiles
from rasterio.merge import merge
from pyproj import Transformer
from sklearn.metrics import silhouette_samples, silhouette_score, ConfusionMatrixDisplay
from rasterstats import zonal_stats
from geocube.api.core import make_geocube
from glob import glob
from sklearn.cluster import KMeans, MiniBatchKMeans
warnings.filterwarnings('ignore')

Code

#ee.Authenticate()

Code

ee.Initialize()

Code

######################## subset L8 collection   ########################################
##################################################################################
def getFactorImg(factorNames, image):
    factorList = image.toDictionary().select(factorNames).values()
    return ee.Image.constant(factorList)

######################## pandas normalise   ########################################
##################################################################################
def NormalizeData(data):
    return (data - np.min(data)) / (np.max(data) - np.min(data))

######################## L8 cloud mask   ########################################
##################################################################################
def getQABits(image, start, end, mask):
    pattern = 0
    for i in range(start,end-1):
        pattern += 2**i
    return image.select([0], [mask]).bitwiseAnd(pattern).rightShift(start)


def maskQuality(image):
    QA = image.select('QA_PIXEL')
    shadowMask = getQABits(QA,3,3,'cloud_shadow')
    cloudMask = getQABits(QA,5,5,'cloud')
    cirrusMask = getQABits(QA,9,9,'cirrus_detected')
    scaleImg = getFactorImg(['REFLECTANCE_MULT_BAND_.|TEMPERATURE_MULT_BAND_ST_B10'], image)
    offsetImg = getFactorImg(['REFLECTANCE_ADD_BAND_.|TEMPERATURE_ADD_BAND_ST_B10'], image)
    scaled = image.select('SR_B.|ST_B10').multiply(scaleImg).add(offsetImg)
    
    return (image.addBands(scaled, None, True)
            .updateMask(cirrusMask)
            # .updateMask(cloudMask)
            # .updateMask(shadowMask)
            )



def prepSrL8(image):
    dilateMask = image.select('QA_PIXEL').bitwiseAnd(int('11111', 0)).eq(0)
    cirrusMask = image.select('QA_PIXEL').bitwiseAnd(int('11111', 2)).eq(0)
    cloudMask = image.select('QA_PIXEL').bitwiseAnd(int('11111', 3)).eq(0)
    shadowMask =  image.select('QA_PIXEL').bitwiseAnd(int('11111', 4)).eq(0)
    snowMask = image.select('QA_PIXEL').bitwiseAnd(int('11111', 5)).eq(0)
    saturationMask = image.select('QA_RADSAT').eq(0)
    scaleImg = getFactorImg(['REFLECTANCE_MULT_BAND_.|TEMPERATURE_MULT_BAND_ST_B10'], image)
    offsetImg = getFactorImg(['REFLECTANCE_ADD_BAND_.|TEMPERATURE_ADD_BAND_ST_B10'], image)
    scaled = image.select('SR_B.|ST_B10').multiply(scaleImg).add(offsetImg)

    return (image.addBands(scaled, None, True)
                 .updateMask(cirrusMask)
                 # .updateMask(dilateMask)
                 # .updateMask(cloudMask)
                 # .updateMask(shadowMask)
                 # .updateMask(snowMask)
                 #.updateMask(saturationMask)
           )

######################## L8 scale factors   ####################################
##################################################################################
def apply_scale_factors(image):
    opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2)
    thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0)
    return image.addBands(opticalBands, None, True).addBands(thermalBands, None, True)

######################## convert L8 to SST   ####################################
##################################################################################
def sst(img):
    thermalBands = img.select('ST_TRAD')
    ln = k2_im.divide((k1_im.divide(thermalBands).add(1).log())).multiply(0.001).clip(site).rename("SST")
    #.convolve(kern)
    return img.addBands(ln)

######################## scale image   ####################################
##################################################################################
def scale(image):
    return ee.Image(image).multiply(0.001)

def addNDAVI(image):
    ndvi = image.normalizedDifference(['SR_B1', 'SR_B5']).rename('NDAVI')
    return image.addBands(ndvi)

def addNDWI(image):
    ndvi = image.normalizedDifference(['SR_B5', 'SR_B6']).rename('NDWI')
    return image.addBands(ndvi)


# Function to add raster to map
def addras(Map, layer, mi, ma, name):
    Map.addLayer(layer,
                 {'bands': ['SR_B4', 'SR_B3', 'SR_B2'],'max': ma, 'min' : mi},
                 name=name)
    
# Function to add raster to map
def addl5(Map, layer, mi, ma, name):
    Map.addLayer(layer,
                 {'bands': ['SR_B3', 'SR_B2', 'SR_B1'],'max': ma, 'min' : mi},
                 name=name)
    
def fmask(image):
    # see https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC09_C02_T1_L2
    # Bit 0 - Fill
    # Bit 1 - Dilated Cloud
    # Bit 2 - Cirrus
    # Bit 3 - Cloud
    # Bit 4 - Cloud Shadow
    qaMask = image.select("QA_PIXEL").bitwiseAnd(int("11111", 2)).eq(0)

    # Apply the scaling factors to the appropriate bands.
    opticalBands = image.select("SR_B.").multiply(0.0000275).add(-0.2)

    # Replace the original bands with the scaled ones and apply the masks.
    return image.addBands(opticalBands, None, True).updateMask(qaMask)

# def addkNDAVI(image):
#     #Compute D2 a rename it to d2
#     var D2 = nir.subtract(red).pow(2)
#     .select([0],['d2']);
#     // Gamma, defined as 1/sigmaˆ2
#     var gamma = ee.Number(4e6).multiply(-2.0);
#     // Compute kernel (k) and KNDVI
#     var k = D2.divide(gamma).exp();
#     var kndvi = ee.Image.constant(1)
#     .subtract(k).divide(
#     ee.Image.constant(1).add(k))
#     .select([0],['knd']);
#     return image.addBands(kndvi);

def landsat5_to_landsat8(landsat5_img):
    #, landsat8_img
    """
    Corrects a Landsat 5 image to Landsat 8 by performing cross-calibration using Earth Engine.
    
    Parameters:
    landsat5_img (ee.Image): The Landsat 5 image to be corrected.
    landsat8_img (ee.Image): The Landsat 8 image to use for cross-calibration.
    
    Returns:
    ee.Image: The corrected Landsat 5 image.
    """
    
    # Define Landsat 5 and Landsat 8 spectral response functions
    l5_sr_func = [0.7661, 1.4475, 1.0049, 0.3734, 0.3362, 0.1446] # Landsat 5
    l8_sr_func = [0.3029, 0.2786, 0.4733, 0.5599, 0.508, 0.1872] # Landsat 8
    
    # Compute conversion factors for each band
    cf = [l8_sr_func[i] / l5_sr_func[i] for i in range(len(l5_sr_func))]
    
    # Apply conversion factors to the Landsat 5 image
    corrected_img = landsat5_img.multiply(cf).toFloat()
    
    # Set the corrected image's metadata to match the Landsat 8 image
    # corrected_img = corrected_img.set('system:time_start', landsat8_img.get('system:time_start'))
    # corrected_img = corrected_img.set('system:time_end', landsat8_img.get('system:time_end'))
    # corrected_img = corrected_img.set('system:footprint', landsat8_img.get('system:footprint'))
    # corrected_img = corrected_img.set('sensor_id', 'OLI_TIRS')
    
    return corrected_img


def export_images_to_cloud(image_list, bucket_name, name_list):
    """
    Exports a list of images to Google Cloud Storage.
    
    Parameters:
    image_list (list): A list of ee.Image objects to export.
    bucket_name (str): The name of the Google Cloud Storage bucket to export the images to.
    
    Returns:
    None
    """
    
    # Loop through the list of images and export each one
    for i, img in enumerate(image_list):
        # Define the output file name
        output_name = name_list[i]
        
        # Export the image to Google Cloud Storage
        task = ee.batch.Export.image.toCloudStorage(image=img,
                                                    description=output_name,
                                                    bucket=bucket_name,
                                                    fileNamePrefix=output_name,
                                                    scale=30,
                                                    crs = 'EPSG:4326',
                                                    fileFormat = 'GeoTIFF',
                                                    maxPixels=1e13)
        
        # Start the export task
        task.start()

Code

sites = geemap.shp_to_ee(os.path.join(os.getcwd(),
          'Data\\ICoAST sites\\WGS84\\ICOAST_midwest.shp'))

kal = geemap.shp_to_ee(os.path.join(os.getcwd(),
          'Data\\ICoAST sites\\Sites\\Kalbarri.shp'))

Code

start_date = '2012-01-01'
end_date = '2012-12-30'
optical_bands = ['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B7']
# L5_2012 = (
#     ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
#     .filterBounds(sites)
#     .filterMetadata('CLOUD_COVER', 'less_than', 5)
#     .filterDate(start_date, end_date)
#     .map(fmask)
#     .median()
#     .clip(sites)

# )


# L5_2012_cor = landsat5_to_landsat8(L5_2012)

start_date = '2011-01-01'
end_date = '2011-12-30'

L5_2011 = (
    ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
    #map(landsat5_to_landsat8)
    .filterBounds(sites)
    .filterMetadata('CLOUD_COVER', 'less_than', 5)
    .filterDate(start_date, end_date)
    .map(fmask)
    .median()
    .clip(sites)
    .select(optical_bands)
)

L5_2011_cor = landsat5_to_landsat8(L5_2011)

start_date = '2010-01-01'
end_date = '2010-12-30'

L5_2010 = (
    ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
    .filterBounds(sites)
    .filterMetadata('CLOUD_COVER', 'less_than', 5)
    .filterDate(start_date, end_date)
    .map(fmask)
    .median()
    .clip(sites)
    .select(optical_bands)
)

L5_2010_cor = landsat5_to_landsat8(L5_2010)

start_date = '2009-01-01'
end_date = '2009-12-30'

L5_2009 = (
    ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
    .filterBounds(sites)
    .filterMetadata('CLOUD_COVER', 'less_than', 5)
    .filterDate(start_date, end_date)
    .map(fmask)
    .median()
    .clip(sites)
    .select(optical_bands)
)

L5_2009_cor = landsat5_to_landsat8(L5_2009)

start_date = '2008-01-01'
end_date = '2008-12-30'

L5_2008 = (
    ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
    .filterBounds(sites)
    .filterMetadata('CLOUD_COVER', 'less_than', 5)
    .filterDate(start_date, end_date)
    .map(fmask)
    .median()
    .clip(sites)
    .select(optical_bands)
)

L5_2008_cor = landsat5_to_landsat8(L5_2008)

start_date = '2007-01-01'
end_date = '2007-12-30'

L5_2007 = (
    ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
    .filterBounds(sites)
    .filterMetadata('CLOUD_COVER', 'less_than', 5)
    .filterDate(start_date, end_date)
    .map(fmask)
    .median()
    .clip(sites)
    .select(optical_bands)
)

L5_2007_cor = landsat5_to_landsat8(L5_2007)


start_date = '2006-01-01'
end_date = '2006-12-30'

L5_2006 = (
    ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
    .filterBounds(sites)
    .filterMetadata('CLOUD_COVER', 'less_than', 5)
    .filterDate(start_date, end_date)
    .map(fmask)
    .median()
    .clip(sites)
    .select(optical_bands)
)

L5_2006_cor = landsat5_to_landsat8(L5_2006)

start_date = '2005-01-01'
end_date = '2005-12-30'

L5_2005 = (
    ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
    .filterBounds(sites)
    .filterMetadata('CLOUD_COVER', 'less_than', 5)
    .filterDate(start_date, end_date)
    .map(fmask)
    .median()
    .clip(sites)
    .select(optical_bands)
)

L5_2005_cor = landsat5_to_landsat8(L5_2005)

start_date = '2004-01-01'
end_date = '2004-12-30'

L5_2004 = (
    ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
    .filterBounds(sites)
    .filterMetadata('CLOUD_COVER', 'less_than', 5)
    .filterDate(start_date, end_date)
    .map(fmask)
    .median()
    .clip(sites)
    .select(optical_bands)
)

L5_2004_cor = landsat5_to_landsat8(L5_2004)

L5 bands Corrected to L8 values

Code

Map = geemap.Map(toolbar_ctrl=True, layer_ctrl=True, add_google_map = False)
#legend = LegendControl({}, name="Landsat 8 2013-2022", position="topleft")
#Map.add_control(legend)
#addl5(Map, L5_2012_cor, 0, 0.05, '2012')
addl5(Map, L5_2011_cor, 0, 0.02, '2011')
addl5(Map, L5_2010_cor, 0, 0.02, '2010')
addl5(Map, L5_2009_cor, 0, 0.02, '2009')
addl5(Map, L5_2008_cor, 0, 0.02, '2008')
addl5(Map, L5_2007_cor, 0, 0.02, '2007')
addl5(Map, L5_2006_cor, 0, 0.02, '2006')
addl5(Map, L5_2005_cor, 0, 0.02, '2005')
addl5(Map, L5_2004_cor, 0, 0.02, '2004')
Map

L5 uncorrected

Code

Map2 = geemap.Map(toolbar_ctrl=True, layer_ctrl=True, add_google_map = False)
#legend = LegendControl({}, name="Landsat 8 2013-2022", position="topleft")
#Map.add_control(legend)
#addl5(Map, L5_2012_cor, 0, 0.05, '2012')
addl5(Map2, L5_2011, 0, 0.05, '2011')
addl5(Map2, L5_2010, 0, 0.05, '2010')
addl5(Map2, L5_2009, 0, 0.05, '2009')
addl5(Map2, L5_2008, 0, 0.05, '2008')
addl5(Map2, L5_2007, 0, 0.05, '2007')
addl5(Map2, L5_2006, 0, 0.05, '2006')
addl5(Map2, L5_2005, 0, 0.05, '2005')
addl5(Map2, L5_2004, 0, 0.05, '2004')
Map2

Code

Map.centerObject(sites.geometry(), 5)
Map

Code

im_list = [L5_2011_cor, L5_2010_cor, L5_2009_cor, L5_2008_cor,
          L5_2007_cor, L5_2006_cor, L5_2005_cor,L5_2004_cor]

im_names = ['L5_2011_cor', 'L5_2010_cor', 'L5_2009_cor', 'L5_2008_cor',
           'L5_2007_cor', 'L5_2006_cor', 'L5_2005_cor','L5_2004_cor']

Code

#export_images_to_cloud(image_list = im_list, bucket_name = 'cog-bucket', name_list = im_names)